This notebook reads in results for the MultiBLUP model using 5000 random gene group to do the following:
- summarize the number of models that pass filtering criteria
- examine the distribution of the likelihood ratio (see Edwards et al 2016)
- establish 95 percentiles for heritability and the likelihood ratio
knitr::opts_knit$set(root.dir = "~/Documents/AA-GenomicPrediction/",
warning = FALSE)
knitr::opts_chunk$set(fig.width = 10)
library(tidyverse)
library(knitr)
library(kableExtra)
library(quantreg)
library(data.table)
library(viridis)
library(emdbook)
library(gridExtra)
library(cowplot)
Size distribution of null pathways
Number of markers should have a uniform distribution
aa_fam <- read.csv("data/processed/aa_family_ids.csv",
header = TRUE)
null_dist <- read_csv("reports/lr_null_results.csv") %>%
# set amino acid family ids
left_join(., aa_fam, by = "trait")
null_dist %>%
distinct(pathway, size.x) %>%
ggplot(., aes(size.x)) +
geom_histogram(binwidth = 1000) +
xlab("Gene group size (number of SNPs)")

Filtering
Filter out values where the REML model did not converge
Set negative heritability estimates to zero
Export summary to S2 table
S2table <- null_dist %>%
group_by(trait) %>%
filter(Component %in% "Share_K1") %>%
mutate("Number of gene groups" = ifelse(SD < 0 | LR < 0, 0, 1),
"Failed to converge" = ifelse(SD < 0, 1, 0),
"Negative LR" = ifelse(LR < 0, 1, 0)) %>%
select("Number of gene groups", "Failed to converge", "Negative LR") %>%
summarize_all(.funs = "sum")
Adding missing grouping variables: `trait`
S2table %>%
write_csv(., "reports/tableS2_null_summary.csv")
head(S2table)
Apply filters to null distribution
null_dist <- null_dist %>%
# mutate(Share = ifelse(Share < 0, 0, Share),
# Share = ifelse(Share > 1, 1, Share)) %>%
rename(size = "size.x",
n_genes = "n_genes.x") %>%
# remove values where SD of heritability is strange (results in inflated LR)
filter(# SD > 0 & LR >= 0 &
Component %in% "Share_K1") %>%
dplyr::select(family, trait, pathway, size, n_genes, group_size, LR, Component, Share, SD)
head(null_dist)
Likelihood ratio
Distribution
Expect LR to be distributed as \(\chi^2_{df 1}\)~ based on Wilks’ Theorem
null_dist %>%
ggplot(., aes(LR)) +
geom_histogram(bins = 50) +
facet_wrap(~trait, scales = "free", ncol = 6) +
theme_bw()

# ggsave("reports/figures/null_histogram.png", height = 12, width = 8)
QQ plots
Maybe some inflation at smaller group sizes (< 10,000 SNPs)
null_dist %>%
# filter(trait %in% c("ala", "ala_PyrFam", "asp_t", "his", "Total")) %>%
mutate(group_size = fct_relevel(group_size, "[0,10000]")) %>%
ggplot(., aes(sample = LR, colour = group_size)) +
scale_colour_manual(values = viridis(6)[1:5]) +
stat_qq(distribution = qchisq, dparams = list(df = 1.5)) +
facet_wrap( ~ trait, scales = "free", ncol = 6) +
geom_abline(slope = 1, intercept = 0) +
theme_bw() +
theme(legend.position = "top")

# ggsave("reports/figures/FigS5.png", height = 12, width = 8)
LR goodness-of-fit tests
Test if LR follows Wilks’ theorem expectation of \(\chi^2\) distribution
- \(\chi^2\) test
- Kolmogorov-Smirnov test (D)
GoF <- null_dist %>%
select(trait, LR) %>%
group_by(trait) %>%
mutate(
# return distance (D) for Kolmogorov-Smirnov goodness-of-fit test
# df = 1
D1 = ks.test(LR, pchisq, df = 1)$statistic,
# df = mixture of 1 and 2
D1_2 = ks.test(LR, rchibarsq(length(LR), 2, mix = 0.5))$statistic,
# df = 2
D2 = ks.test(LR, pchisq, df = 2)$statistic,
# return p-value for chi-sq goodness of fit test
# df = 1
X2_1 = suppressWarnings(chisq.test(LR, rchisq(length(LR), 1))$p.value),
# df = mixture of 1 and 2
X2_1_2 = suppressWarnings(chisq.test(LR, rchibarsq(length(LR), 2, mix = 0.5))$p.value),
# df = 2
X2_2 = suppressWarnings(chisq.test(LR, rchisq(length(LR), 2))$p.value)) %>%
summarise_all(mean) %>%
select(trait, LR, D1, D1_2, D2, X2_1, X2_1_2, X2_2)
# return table of results
GoF %>%
kable(digits = 4) %>%
kable_styling(bootstrap_options = "striped")
| trait |
LR |
D1 |
D1_2 |
D2 |
X2_1 |
X2_1_2 |
X2_2 |
| ala |
0.7433 |
0.0998 |
0.2528 |
0.4018 |
0.2821 |
0.2821 |
0.2821 |
| ala_PyrFam |
3.6045 |
0.7322 |
0.6888 |
0.6492 |
0.2741 |
0.2741 |
0.2741 |
| ala_t |
1.0870 |
0.0310 |
0.1540 |
0.3021 |
0.2731 |
0.2731 |
0.2731 |
| arg |
1.4788 |
0.2762 |
0.1274 |
0.1341 |
0.2631 |
0.2631 |
0.2631 |
| arg_GluFam |
0.3714 |
0.2746 |
0.4298 |
0.5769 |
0.2997 |
0.2997 |
0.2997 |
| arg_t |
0.5028 |
0.1906 |
0.3290 |
0.4820 |
0.2879 |
0.2879 |
0.2879 |
| asp |
0.6665 |
0.1086 |
0.2638 |
0.4059 |
0.2827 |
0.2827 |
0.2827 |
| asp_AspFam |
1.7296 |
0.2077 |
0.0846 |
0.1587 |
0.2611 |
0.2611 |
0.2611 |
| asp_t |
1.7883 |
0.5250 |
0.3986 |
0.2375 |
0.2666 |
0.2666 |
0.2666 |
| AspFam |
0.6107 |
0.1367 |
0.2860 |
0.4300 |
0.2856 |
0.2856 |
0.2856 |
| AspFam_Asp |
1.6519 |
0.1948 |
0.0866 |
0.1760 |
0.2620 |
0.2620 |
0.2620 |
| BCAA |
0.5036 |
0.1833 |
0.3398 |
0.4773 |
0.2855 |
0.2855 |
0.2855 |
| gln |
0.4680 |
0.2112 |
0.3460 |
0.5012 |
0.2872 |
0.2872 |
0.2872 |
| gln_GluFam |
0.3530 |
0.3207 |
0.4626 |
0.6198 |
0.3076 |
0.3076 |
0.3076 |
| gln_t |
0.3872 |
0.2956 |
0.4464 |
0.5948 |
0.3045 |
0.3045 |
0.3045 |
| glu |
0.4454 |
0.2118 |
0.3526 |
0.5077 |
0.2893 |
0.2893 |
0.2893 |
| glu_GluFam |
0.4530 |
0.2315 |
0.3810 |
0.5324 |
0.2956 |
0.2956 |
0.2956 |
| glu_t |
0.9756 |
0.0514 |
0.1918 |
0.3478 |
0.2776 |
0.2776 |
0.2776 |
| GluFam |
0.7638 |
0.1049 |
0.2222 |
0.3325 |
0.2718 |
0.2718 |
0.2718 |
| GluFam_glu |
0.4788 |
0.2183 |
0.3710 |
0.5199 |
0.2935 |
0.2935 |
0.2935 |
| gly |
0.4790 |
0.2193 |
0.3564 |
0.5036 |
0.2823 |
0.2823 |
0.2823 |
| gly_SerFam |
0.3489 |
0.2950 |
0.4408 |
0.5970 |
0.3020 |
0.3020 |
0.3020 |
| gly_t |
0.2541 |
0.3801 |
0.5304 |
0.6819 |
0.3093 |
0.3093 |
0.3093 |
| his |
3.2382 |
0.5027 |
0.4780 |
0.4588 |
0.2677 |
0.2677 |
0.2677 |
| his_GluFam |
1.0107 |
0.0314 |
0.1876 |
0.3334 |
0.2752 |
0.2752 |
0.2752 |
| his_t |
0.7742 |
0.0967 |
0.2296 |
0.3924 |
0.2804 |
0.2804 |
0.2804 |
| ile |
0.5745 |
0.1624 |
0.3138 |
0.4646 |
0.2866 |
0.2866 |
0.2866 |
| ile_AspFam |
0.6032 |
0.1544 |
0.3024 |
0.4568 |
0.2872 |
0.2872 |
0.2872 |
| ile_BCAA |
3.1959 |
0.6695 |
0.6028 |
0.5313 |
0.2671 |
0.2671 |
0.2671 |
| ile_t |
0.8015 |
0.0837 |
0.2388 |
0.3824 |
0.2814 |
0.2814 |
0.2814 |
| leu |
0.5778 |
0.1496 |
0.2926 |
0.4426 |
0.2843 |
0.2843 |
0.2843 |
| leu_BCAA |
0.7253 |
0.0989 |
0.2604 |
0.3997 |
0.2830 |
0.2830 |
0.2830 |
| leu_PyrFam |
1.1848 |
0.1292 |
0.1300 |
0.2752 |
0.2714 |
0.2714 |
0.2714 |
| leu_t |
0.9978 |
0.0251 |
0.1642 |
0.3169 |
0.2748 |
0.2748 |
0.2748 |
| lys |
0.6401 |
0.1277 |
0.2592 |
0.3825 |
0.2787 |
0.2787 |
0.2787 |
| lys_AspFam |
3.0090 |
0.6554 |
0.5930 |
0.5343 |
0.2718 |
0.2718 |
0.2718 |
| lys_t |
0.4535 |
0.2243 |
0.3788 |
0.5265 |
0.2934 |
0.2934 |
0.2934 |
| met |
0.5487 |
0.1559 |
0.2996 |
0.4417 |
0.2834 |
0.2834 |
0.2834 |
| met_AspFam |
-6.3774 |
0.9996 |
0.9996 |
0.9996 |
0.2861 |
0.2861 |
0.2861 |
| met_t |
0.5298 |
0.1675 |
0.3168 |
0.4633 |
0.2865 |
0.2865 |
0.2865 |
| phe |
0.9022 |
0.0539 |
0.2036 |
0.3560 |
0.2788 |
0.2788 |
0.2788 |
| phe_ShikFam |
0.4992 |
0.1990 |
0.3444 |
0.4984 |
0.2903 |
0.2903 |
0.2903 |
| phe_t |
1.0603 |
0.0216 |
0.1622 |
0.3115 |
0.2730 |
0.2730 |
0.2730 |
| pro |
0.4705 |
0.1964 |
0.3382 |
0.4860 |
0.2889 |
0.2889 |
0.2889 |
| pro_GluFam |
0.5255 |
0.1719 |
0.2962 |
0.4496 |
0.2829 |
0.2829 |
0.2829 |
| pro_t |
0.4862 |
0.1919 |
0.3382 |
0.4794 |
0.2871 |
0.2871 |
0.2871 |
| PyrFam |
0.5045 |
0.1926 |
0.3556 |
0.4928 |
0.2909 |
0.2909 |
0.2909 |
| ser |
0.4128 |
0.2313 |
0.3884 |
0.5297 |
0.2930 |
0.2930 |
0.2930 |
| ser_SerFam |
0.4666 |
0.2395 |
0.3872 |
0.5395 |
0.2989 |
0.2989 |
0.2989 |
| ser_t |
0.4228 |
0.2611 |
0.4080 |
0.5560 |
0.3015 |
0.3015 |
0.3015 |
| SerFam |
0.5007 |
0.2011 |
0.3414 |
0.4788 |
0.2827 |
0.2827 |
0.2827 |
| ShikFam |
1.6064 |
0.1484 |
0.0436 |
0.1879 |
0.2628 |
0.2628 |
0.2628 |
| thr |
1.4055 |
0.1325 |
0.0666 |
0.2166 |
0.2657 |
0.2657 |
0.2657 |
| thr_AspFam |
0.7956 |
0.0949 |
0.2530 |
0.3914 |
0.2817 |
0.2817 |
0.2817 |
| total |
0.5574 |
0.1456 |
0.2850 |
0.4249 |
0.2818 |
0.2818 |
0.2818 |
| trp |
0.6488 |
0.1471 |
0.2976 |
0.4487 |
0.2861 |
0.2861 |
0.2861 |
| trp_ShikFam |
0.4828 |
0.2035 |
0.3354 |
0.4928 |
0.2859 |
0.2859 |
0.2859 |
| trp_t |
0.5313 |
0.1832 |
0.3136 |
0.4585 |
0.2813 |
0.2813 |
0.2813 |
| tyr |
0.5001 |
0.1823 |
0.3240 |
0.4780 |
0.2882 |
0.2882 |
0.2882 |
| tyr_ShikFam |
0.4890 |
0.1888 |
0.3278 |
0.4752 |
0.2862 |
0.2862 |
0.2862 |
| tyr_t |
0.8177 |
0.0762 |
0.2144 |
0.3606 |
0.2811 |
0.2811 |
0.2811 |
| val |
0.6674 |
0.1248 |
0.2408 |
0.3726 |
0.2748 |
0.2748 |
0.2748 |
| val_BCAA |
1.3624 |
0.4419 |
0.2864 |
0.2004 |
0.2746 |
0.2746 |
0.2746 |
| val_PyrFam |
0.6507 |
0.1195 |
0.2688 |
0.4137 |
0.2829 |
0.2829 |
0.2829 |
| val_t |
0.3689 |
0.2647 |
0.4112 |
0.5667 |
0.2959 |
0.2959 |
0.2959 |
# plot results
GoF %>%
select(trait, D1, D1_2, D2) %>%
gather(key = df, value = D, -trait) %>%
mutate(df = recode(df, D1 = "1", D1_2 = "1.5", D2 = "2")) %>%
ggplot(., aes(df, D, group = trait)) +
geom_point() +
geom_line() +
xlab("degrees of freedom") +
ylab("Kolmogorov-Smirnov test statistic (D)") +
facet_wrap(~ trait, scales = "fixed", ncol = 6) +
theme_bw() +
theme(legend.position = "none")

# ggsave("reports/figures/FigS4.png", height = 12, width = 10)
GoF %>%
select(trait, D1, D1_2, D2) %>%
gather(key = df, value = D, -trait) %>%
group_by(trait) %>%
mutate(df = recode(df, D1 = "1", D1_2 = "1.5", D2 = "2")) %>%
slice(which.min(D)) %>%
ungroup() %>%
write_csv("reports/goodness_of_fit.csv") %>%
head()
Empirical thresholds
95 percentile for LR
Fit an additive quantile regression to establish 95 percentile cut off for LR at different pathway sizes (number of markers)
## 95% LR threshold for plotting
lr_fit <- NULL
coefs <- NULL
for (i in unique(null_dist$trait)) {
lr_fit[[i]] <- rqss(as.numeric(LR) ~ qss(size, constraint = "I"),
data = null_dist[null_dist$trait==i,], tau = 0.95)
coefs[[i]] <- data.frame(lr_fit[[i]]$qss$size$xyz) %>%
mutate(lr_95 = X2 + lr_fit[[i]]$coef[1],
size = X1) %>%
dplyr::select(size, lr_95)
}
lr_threshold <- rbindlist(coefs, idcol = "trait")
null_dist <- null_dist %>%
left_join(., lr_threshold, by = c("trait", "size")) %>%
mutate(col = LR >= lr_95)
null_dist %>%
as_tibble() %>%
head()
95 percentile for \(h^2\)
Fit an additive quantile regression to establish 95% quantile cut off for proportion of genomic variance explained (H^2) at different pathway sizes (number of markers)
h2_fit <- NULL
h2_coefs <- NULL
for (i in unique(null_dist$trait)) {
h2_fit[[i]] <- rqss(as.numeric(Share) ~ qss(size, constraint = "I"),
data = null_dist[null_dist$trait==i,], tau = 0.977) #FDR corrected
# data = null_dist[null_dist$trait==i,], tau = 0.95)
h2_coefs[[i]] <- data.frame(h2_fit[[i]]$qss$size$xyz) %>%
mutate(h2_95 = X2 + h2_fit[[i]]$coef[1],
size = X1) %>%
dplyr::select(size, h2_95)
}
h2_threshold <- rbindlist(h2_coefs, idcol = "trait") %>%
as_tibble()
head(h2_threshold)
Plot null distribution
Blue dots are pathways that also pass the LR_95 threshold
Solid line is the 95% percentile for proportion of heritability explained
Dashed line is the infinitesimal expectation (each SNP contributes approximately equal amount of variation)

save results
save(lr_fit, h2_fit, file = "reports/null_distribution.RData")
Error in save(lr_fit, h2_fit, file = "reports/null_distribution.RData") :
object ‘h2_fit’ not found
LS0tCnRpdGxlOiAiQUEtR1A6IHByb2Nlc3MgbnVsbCBkaXN0cmlidXRpb24iCm91dHB1dDogCiAgaHRtbF9ub3RlYm9vazoKICAgIHRvYzogeWVzCiAgICB0b2NfZmxvYXQ6IHllcwogICAgY29kZV9mb2xkaW5nOiBoaWRlCmRhdGU6ICJgciBmb3JtYXQoU3lzLnRpbWUoKSwgJyVkICVCLCAlWScpYCIKZWRpdG9yX29wdGlvbnM6IAogIGNodW5rX291dHB1dF90eXBlOiBpbmxpbmUKLS0tCgpUaGlzIG5vdGVib29rIHJlYWRzIGluIHJlc3VsdHMgZm9yIHRoZSBNdWx0aUJMVVAgbW9kZWwgdXNpbmcgNTAwMCByYW5kb20gZ2VuZSBncm91cCB0byBkbyB0aGUgZm9sbG93aW5nOiAgIAoKKiBzdW1tYXJpemUgdGhlIG51bWJlciBvZiBtb2RlbHMgdGhhdCBwYXNzIGZpbHRlcmluZyBjcml0ZXJpYSAgIAoqIGV4YW1pbmUgdGhlIGRpc3RyaWJ1dGlvbiBvZiB0aGUgbGlrZWxpaG9vZCByYXRpbyAoc2VlIEVkd2FyZHMgZXQgYWwgMjAxNikgICAKKiBlc3RhYmxpc2ggOTUgcGVyY2VudGlsZXMgZm9yIGhlcml0YWJpbGl0eSBhbmQgdGhlIGxpa2VsaWhvb2QgcmF0aW8gIAoKYGBge3Igc2V0dXAsIHJlc3VsdHMgPSAiaGlkZSJ9CmtuaXRyOjpvcHRzX2tuaXQkc2V0KHJvb3QuZGlyID0gIn4vRG9jdW1lbnRzL0FBLUdlbm9taWNQcmVkaWN0aW9uLyIsCiAgICAgICAgICAgICAgICAgICAgIHdhcm5pbmcgPSBGQUxTRSkKa25pdHI6Om9wdHNfY2h1bmskc2V0KGZpZy53aWR0aCA9IDEwKQoKbGlicmFyeSh0aWR5dmVyc2UpCmxpYnJhcnkoa25pdHIpCmxpYnJhcnkoa2FibGVFeHRyYSkKbGlicmFyeShxdWFudHJlZykKbGlicmFyeShkYXRhLnRhYmxlKQpsaWJyYXJ5KHZpcmlkaXMpCmxpYnJhcnkoZW1kYm9vaykKbGlicmFyeShncmlkRXh0cmEpCmxpYnJhcnkoY293cGxvdCkKYGBgCgojIFNpemUgZGlzdHJpYnV0aW9uIG9mIG51bGwgcGF0aHdheXMKTnVtYmVyIG9mIG1hcmtlcnMgc2hvdWxkIGhhdmUgYSB1bmlmb3JtIGRpc3RyaWJ1dGlvbgoKYGBge3IgbG9hZCBkYXRhLCBmaWcud2lkdGggPSA1LCByZXN1bHRzID0gImhpZGUifSAKYWFfZmFtIDwtIHJlYWQuY3N2KCJkYXRhL3Byb2Nlc3NlZC9hYV9mYW1pbHlfaWRzLmNzdiIsCiAgICAgICAgICAgICAgICAgICBoZWFkZXIgPSBUUlVFKQoKbnVsbF9kaXN0IDwtIHJlYWRfY3N2KCJyZXBvcnRzL2xyX251bGxfcmVzdWx0cy5jc3YiKSAlPiUKICAjIHNldCBhbWlubyBhY2lkIGZhbWlseSBpZHMgCiAgbGVmdF9qb2luKC4sIGFhX2ZhbSwgYnkgPSAidHJhaXQiKQoKbnVsbF9kaXN0ICU+JQogIGRpc3RpbmN0KHBhdGh3YXksIHNpemUueCkgJT4lCiAgZ2dwbG90KC4sIGFlcyhzaXplLngpKSArCiAgZ2VvbV9oaXN0b2dyYW0oYmlud2lkdGggPSAxMDAwKSArCiAgeGxhYigiR2VuZSBncm91cCBzaXplIChudW1iZXIgb2YgU05QcykiKQpgYGAKCiMgRmlsdGVyaW5nCkZpbHRlciBvdXQgdmFsdWVzIHdoZXJlIHRoZSBSRU1MIG1vZGVsIGRpZCBub3QgY29udmVyZ2UgIAoKU2V0IG5lZ2F0aXZlIGhlcml0YWJpbGl0eSBlc3RpbWF0ZXMgdG8gemVybyAKCkV4cG9ydCBzdW1tYXJ5IHRvIFMyIHRhYmxlCgpgYGB7ciB0YWJsZSBTMn0KUzJ0YWJsZSA8LSBudWxsX2Rpc3QgJT4lCiAgZ3JvdXBfYnkodHJhaXQpICU+JQogIGZpbHRlcihDb21wb25lbnQgJWluJSAiU2hhcmVfSzEiKSAlPiUKICBtdXRhdGUoIk51bWJlciBvZiBnZW5lIGdyb3VwcyIgPSBpZmVsc2UoU0QgPCAwIHwgTFIgPCAwLCAwLCAxKSwKICAgICAgICAgIkZhaWxlZCB0byBjb252ZXJnZSIgPSBpZmVsc2UoU0QgPCAwLCAxLCAwKSwKICAgICAgICAgIk5lZ2F0aXZlIExSIiA9IGlmZWxzZShMUiA8IDAsIDEsIDApKSAlPiUKICBzZWxlY3QoIk51bWJlciBvZiBnZW5lIGdyb3VwcyIsICJGYWlsZWQgdG8gY29udmVyZ2UiLCAiTmVnYXRpdmUgTFIiKSAlPiUKICBzdW1tYXJpemVfYWxsKC5mdW5zID0gInN1bSIpIAoKUzJ0YWJsZSAlPiUKICB3cml0ZV9jc3YoLiwgInJlcG9ydHMvdGFibGVTMl9udWxsX3N1bW1hcnkuY3N2IikKCmhlYWQoUzJ0YWJsZSkKYGBgCgpBcHBseSBmaWx0ZXJzIHRvIG51bGwgZGlzdHJpYnV0aW9uIAoKYGBge3IgZmlsdGVyIG51bGx9Cm51bGxfZGlzdCA8LSBudWxsX2Rpc3QgJT4lIAogICMgbXV0YXRlKFNoYXJlID0gaWZlbHNlKFNoYXJlIDwgMCwgMCwgU2hhcmUpLAogICMgICAgICAgIFNoYXJlID0gaWZlbHNlKFNoYXJlID4gMSwgMSwgU2hhcmUpKSAlPiUKICByZW5hbWUoc2l6ZSA9ICJzaXplLngiLAogICAgICAgICBuX2dlbmVzID0gIm5fZ2VuZXMueCIpICU+JSAKICAjIHJlbW92ZSB2YWx1ZXMgd2hlcmUgU0Qgb2YgaGVyaXRhYmlsaXR5IGlzIHN0cmFuZ2UgKHJlc3VsdHMgaW4gaW5mbGF0ZWQgTFIpCiAgZmlsdGVyKCMgU0QgPiAwICYgTFIgPj0gMCAmCiAgICBDb21wb25lbnQgJWluJSAiU2hhcmVfSzEiKSAlPiUKICBkcGx5cjo6c2VsZWN0KGZhbWlseSwgdHJhaXQsIHBhdGh3YXksIHNpemUsIG5fZ2VuZXMsIGdyb3VwX3NpemUsIExSLCBDb21wb25lbnQsIFNoYXJlLCBTRCkgCgpoZWFkKG51bGxfZGlzdCkKYGBgCgojIExpa2VsaWhvb2QgcmF0aW8KCiMjIERpc3RyaWJ1dGlvbgpFeHBlY3QgTFIgdG8gYmUgZGlzdHJpYnV0ZWQgYXMgJFxjaGleMl97ZGYgMX0kfiBiYXNlZCBvbiBXaWxrcycgVGhlb3JlbSAKCmBgYHtyIG51bGwgbHIgZGlzdHJpYnV0aW9uLCBmaWcuaGVpZ2h0ID0gMTQsIGZpZy53aWR0aCA9IDh9Cm51bGxfZGlzdCAlPiUKICAgIGdncGxvdCguLCBhZXMoTFIpKSArIAogICAgZ2VvbV9oaXN0b2dyYW0oYmlucyA9IDUwKSArCiAgZmFjZXRfd3JhcCh+dHJhaXQsIHNjYWxlcyA9ICJmcmVlIiwgbmNvbCA9IDYpICsKICB0aGVtZV9idygpIAoKIyBnZ3NhdmUoInJlcG9ydHMvZmlndXJlcy9udWxsX2hpc3RvZ3JhbS5wbmciLCBoZWlnaHQgPSAxMiwgd2lkdGggPSA4KQpgYGAKCiMjIFFRIHBsb3RzCgpNYXliZSBzb21lIGluZmxhdGlvbiBhdCBzbWFsbGVyIGdyb3VwIHNpemVzICg8IDEwLDAwMCBTTlBzKQoKYGBge3IgcXFjaGlzcSwgZmlnLmhlaWdodCA9IDE0LCBmaWcud2lkdGggPSA4fQpudWxsX2Rpc3QgJT4lCiAgIyBmaWx0ZXIodHJhaXQgJWluJSBjKCJhbGEiLCAiYWxhX1B5ckZhbSIsICJhc3BfdCIsICJoaXMiLCAiVG90YWwiKSkgJT4lIAogIG11dGF0ZShncm91cF9zaXplID0gZmN0X3JlbGV2ZWwoZ3JvdXBfc2l6ZSwgIlswLDEwMDAwXSIpKSAlPiUKICBnZ3Bsb3QoLiwgYWVzKHNhbXBsZSA9IExSLCBjb2xvdXIgPSBncm91cF9zaXplKSkgKyAKICBzY2FsZV9jb2xvdXJfbWFudWFsKHZhbHVlcyA9IHZpcmlkaXMoNilbMTo1XSkgKyAKICBzdGF0X3FxKGRpc3RyaWJ1dGlvbiA9IHFjaGlzcSwgZHBhcmFtcyA9IGxpc3QoZGYgPSAxLjUpKSArCiAgZmFjZXRfd3JhcCggfiB0cmFpdCwgc2NhbGVzID0gImZyZWUiLCBuY29sID0gNikgKwogIGdlb21fYWJsaW5lKHNsb3BlID0gMSwgaW50ZXJjZXB0ID0gMCkgKwogIHRoZW1lX2J3KCkgKyAKICB0aGVtZShsZWdlbmQucG9zaXRpb24gPSAidG9wIikgCgojIGdnc2F2ZSgicmVwb3J0cy9maWd1cmVzL0ZpZ1M1LnBuZyIsIGhlaWdodCA9IDEyLCB3aWR0aCA9IDgpCmBgYAoKIyMgTFIgZ29vZG5lc3Mtb2YtZml0IHRlc3RzIApUZXN0IGlmIExSIGZvbGxvd3MgV2lsa3MnIHRoZW9yZW0gZXhwZWN0YXRpb24gb2YgJFxjaGleMiQgZGlzdHJpYnV0aW9uICAgCgoqICRcY2hpXjIkIHRlc3QgICAgCiogS29sbW9nb3Jvdi1TbWlybm92IHRlc3QgKEQpICAKCmBgYHtyIHRlc3QgZ29vZG5lc3Mgb2YgZml0LCByZXN1bHRzID0gImhpZGUiLCBtZXNzYWdlID0gRkFMU0V9CkdvRiA8LSBudWxsX2Rpc3QgJT4lIAogIHNlbGVjdCh0cmFpdCwgTFIpICU+JQogIGdyb3VwX2J5KHRyYWl0KSAlPiUKICBtdXRhdGUoCiAgICAjIHJldHVybiBkaXN0YW5jZSAoRCkgZm9yIEtvbG1vZ29yb3YtU21pcm5vdiBnb29kbmVzcy1vZi1maXQgdGVzdCAKICAgICMgZGYgPSAxCiAgICBEMSA9IGtzLnRlc3QoTFIsIHBjaGlzcSwgZGYgPSAxKSRzdGF0aXN0aWMsCiAgICAjIGRmID0gbWl4dHVyZSBvZiAxIGFuZCAyIAogICAgRDFfMiA9IGtzLnRlc3QoTFIsIHJjaGliYXJzcShsZW5ndGgoTFIpLCAyLCBtaXggPSAwLjUpKSRzdGF0aXN0aWMsCiAgICAjIGRmID0gMgogICAgRDIgPSBrcy50ZXN0KExSLCBwY2hpc3EsIGRmID0gMikkc3RhdGlzdGljLAogICAgIyByZXR1cm4gcC12YWx1ZSBmb3IgY2hpLXNxIGdvb2RuZXNzIG9mIGZpdCB0ZXN0IAogICAgIyBkZiA9IDEKICAgIFgyXzEgPSBzdXBwcmVzc1dhcm5pbmdzKGNoaXNxLnRlc3QoTFIsIHJjaGlzcShsZW5ndGgoTFIpLCAxKSkkcC52YWx1ZSksCiAgICAjIGRmID0gbWl4dHVyZSBvZiAxIGFuZCAyCiAgICBYMl8xXzIgPSBzdXBwcmVzc1dhcm5pbmdzKGNoaXNxLnRlc3QoTFIsIHJjaGliYXJzcShsZW5ndGgoTFIpLCAyLCBtaXggPSAwLjUpKSRwLnZhbHVlKSwKICAgICMgZGYgPSAyCiAgICBYMl8yID0gc3VwcHJlc3NXYXJuaW5ncyhjaGlzcS50ZXN0KExSLCByY2hpc3EobGVuZ3RoKExSKSwgMikpJHAudmFsdWUpKSAlPiUKICBzdW1tYXJpc2VfYWxsKG1lYW4pICU+JQogIHNlbGVjdCh0cmFpdCwgTFIsIEQxLCBEMV8yLCBEMiwgWDJfMSwgWDJfMV8yLCBYMl8yKQpgYGAKCmBgYHtyIGdvZiB0YWJsZX0KIyByZXR1cm4gdGFibGUgb2YgcmVzdWx0cyAKR29GICU+JSAKICBrYWJsZShkaWdpdHMgPSA0KSAlPiUKICBrYWJsZV9zdHlsaW5nKGJvb3RzdHJhcF9vcHRpb25zID0gInN0cmlwZWQiKQpgYGAKCmBgYHtyIHBsb3QgS1Mgc3RhdGlzdGljLCBmaWcuaGVpZ2h0ID0gMTQsIGZpZy53aWR0aCA9IDh9CiMgcGxvdCByZXN1bHRzCkdvRiAlPiUKICBzZWxlY3QodHJhaXQsIEQxLCBEMV8yLCBEMikgJT4lCiAgZ2F0aGVyKGtleSA9IGRmLCB2YWx1ZSA9IEQsIC10cmFpdCkgJT4lCiAgbXV0YXRlKGRmID0gcmVjb2RlKGRmLCBEMSA9ICIxIiwgRDFfMiA9ICIxLjUiLCBEMiA9ICIyIikpICU+JQogIGdncGxvdCguLCBhZXMoZGYsIEQsIGdyb3VwID0gdHJhaXQpKSArCiAgZ2VvbV9wb2ludCgpICsgCiAgZ2VvbV9saW5lKCkgKwogIHhsYWIoImRlZ3JlZXMgb2YgZnJlZWRvbSIpICsKICB5bGFiKCJLb2xtb2dvcm92LVNtaXJub3YgdGVzdCBzdGF0aXN0aWMgKEQpIikgKyAKICBmYWNldF93cmFwKH4gdHJhaXQsIHNjYWxlcyA9ICJmaXhlZCIsIG5jb2wgPSA2KSArCiAgdGhlbWVfYncoKSArCiAgdGhlbWUobGVnZW5kLnBvc2l0aW9uID0gIm5vbmUiKQoKIyBnZ3NhdmUoInJlcG9ydHMvZmlndXJlcy9GaWdTNC5wbmciLCBoZWlnaHQgPSAxMiwgd2lkdGggPSAxMCkgIApgYGAKCmBgYHtyIGdldCBtaW5pbXVtIEQgdmFsdWV9CkdvRiAlPiUKICBzZWxlY3QodHJhaXQsIEQxLCBEMV8yLCBEMikgJT4lCiAgZ2F0aGVyKGtleSA9IGRmLCB2YWx1ZSA9IEQsIC10cmFpdCkgJT4lCiAgZ3JvdXBfYnkodHJhaXQpICU+JQogIG11dGF0ZShkZiA9IHJlY29kZShkZiwgRDEgPSAiMSIsIEQxXzIgPSAiMS41IiwgRDIgPSAiMiIpKSAlPiUKICBzbGljZSh3aGljaC5taW4oRCkpICU+JQogIHVuZ3JvdXAoKSAlPiUKICB3cml0ZV9jc3YoInJlcG9ydHMvZ29vZG5lc3Nfb2ZfZml0LmNzdiIpICU+JQogIGhlYWQoKQpgYGAKCiMgRW1waXJpY2FsIHRocmVzaG9sZHMKCiMjIDk1IHBlcmNlbnRpbGUgZm9yIExSIAoKRml0IGFuIGFkZGl0aXZlIHF1YW50aWxlIHJlZ3Jlc3Npb24gdG8gZXN0YWJsaXNoIDk1IHBlcmNlbnRpbGUgY3V0IG9mZiBmb3IgTFIgYXQgZGlmZmVyZW50IHBhdGh3YXkgc2l6ZXMgKG51bWJlciBvZiBtYXJrZXJzKSAKCmBgYHtyIDk1IHF1YW50aWxlIGZvciBMUn0KIyMgOTUlIExSIHRocmVzaG9sZCBmb3IgcGxvdHRpbmcKbHJfZml0IDwtIE5VTEwKY29lZnMgPC0gTlVMTAoKZm9yIChpIGluIHVuaXF1ZShudWxsX2Rpc3QkdHJhaXQpKSB7CiAgbHJfZml0W1tpXV0gPC0gcnFzcyhhcy5udW1lcmljKExSKSB+IHFzcyhzaXplLCBjb25zdHJhaW50ID0gIkkiKSwgCiAgICAgICAgICAgICAgICAgICAgICBkYXRhID0gbnVsbF9kaXN0W251bGxfZGlzdCR0cmFpdD09aSxdLCB0YXUgPSAwLjk1KSAKICBjb2Vmc1tbaV1dIDwtIGRhdGEuZnJhbWUobHJfZml0W1tpXV0kcXNzJHNpemUkeHl6KSAlPiUKICAgIG11dGF0ZShscl85NSA9IFgyICsgbHJfZml0W1tpXV0kY29lZlsxXSwKICAgICAgICAgICBzaXplID0gWDEpICU+JQogICAgZHBseXI6OnNlbGVjdChzaXplLCBscl85NSkKfQoKbHJfdGhyZXNob2xkIDwtIHJiaW5kbGlzdChjb2VmcywgaWRjb2wgPSAidHJhaXQiKQoKbnVsbF9kaXN0IDwtIG51bGxfZGlzdCAlPiUKICBsZWZ0X2pvaW4oLiwgbHJfdGhyZXNob2xkLCBieSA9IGMoInRyYWl0IiwgInNpemUiKSkgJT4lCiAgbXV0YXRlKGNvbCA9IExSID49IGxyXzk1KSAKCm51bGxfZGlzdCAlPiUKICBhc190aWJibGUoKSAlPiUKICBoZWFkKCkKYGBgCgojIyA5NSBwZXJjZW50aWxlIGZvciAkaF4yJCAKCkZpdCBhbiBhZGRpdGl2ZSBxdWFudGlsZSByZWdyZXNzaW9uIHRvIGVzdGFibGlzaCA5NSUgcXVhbnRpbGUgY3V0IG9mZiBmb3IgcHJvcG9ydGlvbiBvZiBnZW5vbWljIHZhcmlhbmNlIGV4cGxhaW5lZCAoSF4yKSBhdCBkaWZmZXJlbnQgcGF0aHdheSBzaXplcyAobnVtYmVyIG9mIG1hcmtlcnMpIAoKYGBge3IgOTUgcXVhbnRpbGUgZm9yIEgyfQpoMl9maXQgPC0gTlVMTApoMl9jb2VmcyA8LSBOVUxMCgpmb3IgKGkgaW4gdW5pcXVlKG51bGxfZGlzdCR0cmFpdCkpIHsKICBoMl9maXRbW2ldXSA8LSBycXNzKGFzLm51bWVyaWMoU2hhcmUpIH4gcXNzKHNpemUsIGNvbnN0cmFpbnQgPSAiSSIpLCAKICAgICAgICAgICAgICAgICAgICAgIGRhdGEgPSBudWxsX2Rpc3RbbnVsbF9kaXN0JHRyYWl0PT1pLF0sIHRhdSA9IDAuOTc3KSAjRkRSIGNvcnJlY3RlZCAKICAgICAgICAgICAgICAgICAgICAgICMgZGF0YSA9IG51bGxfZGlzdFtudWxsX2Rpc3QkdHJhaXQ9PWksXSwgdGF1ID0gMC45NSkKICBoMl9jb2Vmc1tbaV1dIDwtIGRhdGEuZnJhbWUoaDJfZml0W1tpXV0kcXNzJHNpemUkeHl6KSAlPiUKICAgIG11dGF0ZShoMl85NSA9IFgyICsgaDJfZml0W1tpXV0kY29lZlsxXSwKICAgICAgICAgICBzaXplID0gWDEpICU+JQogICAgZHBseXI6OnNlbGVjdChzaXplLCBoMl85NSkKfQoKaDJfdGhyZXNob2xkIDwtIHJiaW5kbGlzdChoMl9jb2VmcywgaWRjb2wgPSAidHJhaXQiKSAlPiUKICBhc190aWJibGUoKQogIApoZWFkKGgyX3RocmVzaG9sZCkKYGBgCgojIyBQbG90IG51bGwgZGlzdHJpYnV0aW9uCkJsdWUgZG90cyBhcmUgcGF0aHdheXMgdGhhdCBhbHNvIHBhc3MgdGhlIExSXzk1IHRocmVzaG9sZCAgClNvbGlkIGxpbmUgaXMgdGhlIDk1JSBwZXJjZW50aWxlIGZvciBwcm9wb3J0aW9uIG9mIGhlcml0YWJpbGl0eSBleHBsYWluZWQgIApEYXNoZWQgbGluZSBpcyB0aGUgaW5maW5pdGVzaW1hbCBleHBlY3RhdGlvbiAoZWFjaCBTTlAgY29udHJpYnV0ZXMgYXBwcm94aW1hdGVseSBlcXVhbCBhbW91bnQgb2YgdmFyaWF0aW9uKQoKYGBge3IgaDIgOTUgdGhyZXNob2xkLCBmaWcuaGVpZ2h0ID0gMTQsIGZpZy53aWR0aCA9IDh9Cm51bGxfZGlzdCAlPiUKICBnZ3Bsb3QoLiwgYWVzKHNpemUsIFNoYXJlLCBjb2xvdXIgPSBjb2wpKSArCiAgc2NhbGVfY29sb3VyX21hbnVhbCh2YWx1ZXMgPSBjKCJsaWdodGdyZXkiLCAiIzFjOTA5OSIpLAogICAgICAgICAgICAgICAgICAgICAgbGFiZWxzID0gYyhleHByZXNzaW9uKHBhc3RlKExSIDwgTFJbOTVdKSksIAogICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICBleHByZXNzaW9uKHBhc3RlKExSID4gTFJbOTVdKSkpKSArIAogIGdlb21fcG9pbnQocG9zaXRpb24gPSBwb3NpdGlvbl9kb2RnZSh3aWR0aCA9IDAuMyksIGFlcyhzaXplID0gY29sKSwgYWxwaGEgPSAwLjc1KSArIAogIHNjYWxlX3NpemVfbWFudWFsKGd1aWRlID0gIm5vbmUiLCB2YWx1ZXMgPSBjKDAuMSwgMC41KSkgKyAKICBnZW9tX3F1YW50aWxlKG1ldGhvZCA9ICJycXNzIiwgbGFtYmRhID0gNTAwLCBxdWFudGlsZXMgPSBjKDAuOTUpLAogICAgICAgICAgICAgICAgZm9ybXVsYSA9IGFzLmZvcm11bGEoeSB+IHFzcyh4LCBjb25zdHJhaW50ID0gIkkiKSksCiAgICAgICAgICAgICAgICBhZXMobGluZXR5cGUgPSBmYWN0b3IoLi5xdWFudGlsZS4uKSksCiAgICAgICAgICAgICAgICBjb2xvdXIgPSAiYmxhY2siLCBzaG93LmxlZ2VuZCA9IEZBTFNFKSArIAogIGdlb21fYWJsaW5lKGFlcyhpbnRlcmNlcHQgPSAwLCBzbG9wZSA9IDEvMTk5NDUyLCBsaW5ldHlwZSA9ICJkYXNoZWQiKSkgKyAKICBzY2FsZV9saW5ldHlwZV9tYW51YWwobGFiZWxzID0gYyhleHByZXNzaW9uKHBhc3RlKEhbOTVdXjIpKSwgCiAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgImluZmluaXRlc2ltYWwiKSwKICAgICAgICAgICAgICAgICAgICAgICAgdmFsdWVzID0gYygxLDIpKSArCiAgZmFjZXRfd3JhcCggfiB0cmFpdCwgbmNvbCA9IDYpICsKICB4bGltKDAsIDUwMDAwKSArIAogIHRoZW1lX2J3KCkgKwogIHRoZW1lKGF4aXMudGV4dC54ID0gZWxlbWVudF90ZXh0KGFuZ2xlID0gOTAsIGhqdXN0ID0gMSksCiAgICAgICAgbGVnZW5kLnRpdGxlID0gZWxlbWVudF9ibGFuaygpLAogICAgICAgIGxlZ2VuZC50ZXh0LmFsaWduID0gMCkKCiMgZ2dzYXZlKCJyZXBvcnRzL2ZpZ3VyZXMvUzJGaWcudGlmZiIsIGhlaWdodCA9IDEyLCB3aWR0aCA9IDEwKSAgCmBgYAoKIyBzYXZlIHJlc3VsdHMKYGBge3Igc2F2ZSByZXN1bHRzfQpzYXZlKGxyX2ZpdCwgaDJfZml0LCBmaWxlID0gInJlcG9ydHMvbnVsbF9kaXN0cmlidXRpb24uUkRhdGEiKQpgYGAK